*! covtest v0.3 2019-01-03
*! based on Yu, Zhang and  Zuo (2019). Multiple Switching and Data Quality in the Multiple Price List.
capture program drop covtest
program define covtest, rclass sortpreserve byable(recall)
	version 14.2
	syntax varlist(max=2 min=2 numeric) [if] [in], BY(varname) [Level(cilevel)]
	tokenize `varlist'
	args x y treatmentvar
	
	marksample touse
	quietly count if `touse'
	if `r(N)' == 0 {
		error 2000
	}

	* check treatment is binary
	quietly levelsof `by'
	if r(levels) != "0 1" {
		display as error "{bf:`by'} must be binary with levels 0 and 1"
		exit
	}
	

	*** declare temporary variables and macro names
	tempvar x_std y_std w wmssq
	tempname X_0_mean X_0_var Y_0_mean Y_0_var XYcor_0 ///
		X_1_mean X_1_var Y_1_mean Y_1_var XYcor_1 ///
		CorrS_0 CorrS_1 ///
		N_0 N_1 S_0 W_0_mean S_1 W_1_mean VW_0 VW_1 ///
		VS_0 VS_1 P SP P_0 P_1 R_0_fisher R_1_fisher R_fisher ///
		P_CORR_0 P_CORR_1 P_CORR result_matrix
	
	*** do calculations
	quietly summarize `x' if `by'==0 & `touse'
	scalar X_0_mean = r(mean)
	scalar X_0_var = r(Var)
	quietly summarize `y' if `by'==0 & `touse'
	scalar Y_0_mean = r(mean)
	scalar Y_0_var = r(Var)
	quietly corr `x' `y' if `by'==0 & `touse'
	scalar XYcor_0 = r(rho)
	
	quietly summarize `x' if `by'==1 & `touse'
	scalar X_1_mean = r(mean)
	scalar X_1_var = r(Var)
	quietly summarize `y' if `by'==1 & `touse'
	scalar Y_1_mean = r(mean)
	scalar Y_1_var = r(Var)
	quietly corr `x' `y' if `by'==1 & `touse'
	scalar XYcor_1 = r(rho)
	
	quietly generate `x_std' = `x' - X_0_mean if `by'==0 & `touse'
	quietly replace `x_std' = `x' - X_1_mean if `by'==1 & `touse'
	quietly generate `y_std' = `y' - Y_0_mean if `by'==0 & `touse'
	quietly replace `y_std' = `y' - Y_1_mean if `by'==1 & `touse'
	quietly generate `w' = `x_std' * `y_std'
	
	quietly summarize `by' if `by'==0 & `touse'
	scalar N_0 = r(N)
	quietly summarize `by' if `by'==1 & `touse'
	scalar N_1 = r(N)
	quietly summarize `w' if `by'==0 & `touse'
	scalar S_0 = r(sum)/(N_0-1)
	scalar W_0_mean = r(mean)
	quietly summarize `w' if `by'==1 & `touse'
	scalar S_1 = r(sum)/(N_1-1)
	scalar W_1_mean = r(mean)
	
	scalar CorrS_0 = S_0/(sqrt(X_0_var)*sqrt(Y_0_var))
	scalar CorrS_1 = S_1/(sqrt(X_1_var)*sqrt(Y_1_var))
	
	quietly generate `wmssq' = (`w' - S_0)^2 if `by'==0
	quietly replace `wmssq' = (`w' - S_1)^2 if `by'==1
	quietly summarize `wmssq' if `by'==0 & `touse'
	scalar VW_0 = r(mean) - (W_0_mean - S_0)^2
	quietly summarize `wmssq' if `by'==1 & `touse'
	scalar VW_1 = r(mean) - (W_1_mean - S_1)^2
	
	scalar VS_0 = VW_0/N_0 + (X_0_var*Y_0_var*(XYcor_0^2+1))/(N_0*(N_0-1))
	scalar VS_1 = VW_1/N_1 + (X_1_var*Y_1_var*(XYcor_1^2+1))/(N_1*(N_1-1))
	
	local LL_0_95 = S_0 - invnormal(1 - 0.05/2)*sqrt(VS_0)
	local UL_0_95 = S_0 + invnormal(1 - 0.05/2)*sqrt(VS_0)
	local LL_0_99 = S_0 - invnormal(1 - 0.01/2)*sqrt(VS_0)
	local UL_0_99 = S_0 + invnormal(1 - 0.01/2)*sqrt(VS_0)
	scalar P_0 = 2-2*normal(abs(S_0)/sqrt(VS_0))

	local LL_1_95 = S_1 - invnormal(1 - 0.05/2)*sqrt(VS_1)
	local UL_1_95 = S_1 + invnormal(1 - 0.05/2)*sqrt(VS_1)
	local LL_1_99 = S_1 - invnormal(1 - 0.01/2)*sqrt(VS_1)
	local UL_1_99 = S_1 + invnormal(1 - 0.01/2)*sqrt(VS_1)
	scalar P_1 = 2-2*normal(abs(S_1)/sqrt(VS_1))
	
	local LL_95 = (S_0 - S_1) - invnormal(1 - 0.05/2)*sqrt(VS_0 + VS_1)
	local UL_95 = (S_0 - S_1) + invnormal(1 - 0.05/2)*sqrt(VS_0 + VS_1)
	local LL_99 = (S_0 - S_1) - invnormal(1 - 0.01/2)*sqrt(VS_0 + VS_1)
	local UL_99 = (S_0 - S_1) + invnormal(1 - 0.01/2)*sqrt(VS_0 + VS_1)
	scalar P = normal((S_0 - S_1)/sqrt(VS_0 + VS_1))
	scalar SP = 2-2*normal(abs(S_0 - S_1)/sqrt(VS_0 + VS_1))
	
	scalar R_0_fisher = (ln(1+CorrS_0) - ln(1-CorrS_0))/2
	scalar R_1_fisher = (ln(1+CorrS_1) - ln(1-CorrS_1))/2
	scalar R_fisher = R_0_fisher - R_1_fisher
	scalar P_CORR_0 = 2 - 2*normal(abs(R_0_fisher)*sqrt(N_0-3))
	scalar P_CORR_1 = 2 - 2*normal(abs(R_1_fisher)*sqrt(N_1-3))
	scalar P_CORR = 2 - 2*normal(abs(R_fisher)/sqrt(1/(N_0-3) + 1/(N_1-3)))
	
	*** construct result table
	scalar SE_0 = sqrt(VS_0)
	scalar SE_1 = sqrt(VS_1)
	scalar SE = sqrt(VS_0 + VS_1)

	scalar LL_0 = S_0 - invnormal(1 - (1 - 0.01*`level')/2)*sqrt(VS_0)
	scalar UL_0 = S_0 + invnormal(1 - (1 - 0.01*`level')/2)*sqrt(VS_0)
	scalar LL_1 = S_1 - invnormal(1 - (1 - 0.01*`level')/2)*sqrt(VS_1)
	scalar UL_1 = S_1 + invnormal(1 - (1 - 0.01*`level')/2)*sqrt(VS_1)
	scalar LL = (S_0 - S_1) - invnormal(1 - (1 - 0.01*`level')/2)*sqrt(VS_0 + VS_1)
	scalar UL = (S_0 - S_1) + invnormal(1 - (1 - 0.01*`level')/2)*sqrt(VS_0 + VS_1)
	
	local byvar8 = abbrev("`by'",8)
	
	display ""
	display "{text}Test equality of covariance and correlation of {bf:`x'} and {bf:`y'} grouped by {bf:`by'}"
	display "{text}{hline 12}{c TT}{hline 69}{c TT}{hline 23}"
	display "{text}{space 12}{c |}{ralign 6:Obs}  {ralign 11:Covariance}" ///
		" {ralign 11:Std. Err.}{ralign 11:p-value}{ralign 26:[`level'% Conf. Interval]} " ///
		"{c |} {ralign 11:Correlation}{ralign 11:p-value}"
	display "{hline 12}{c +}{hline 69}{c +}{hline 23}"
	display "{text}{ralign 10:0.`byvar8'}" "  {c |}{result}" %6.0g N_0 "     " %8.6g S_0 "    " ///
		%8.6g SE_0 "   " %8.6g P_0 "     " %8.6g LL_0 "     " %8.6g UL_0 " {c |}    " ///
		%8.6g CorrS_0 "   " %8.6g P_CORR_0 "   "
	display "{text}{ralign 10:1.`byvar8'}" "  {c |}{result}" %6.0g N_1 "     " %8.6g S_1 "    " ///
		%8.6g SE_1 "   " %8.6g P_1 "     " %8.6g LL_1 "     " %8.6g UL_1 " {c |}    " ///
		%8.6g CorrS_1 "   " %8.6g P_CORR_1 "   "
	display "{text}{hline 12}{c +}{hline 69}{c +}{hline 23}"
	display "{text}{ralign 10:diff}  {c |}{result}" %6.0g N_0+N_1 "     " %8.6g S_0-S_1 "    " ///
		%8.6g SE "   " %8.6g SP "     " %8.6g LL "     " %8.6g UL " {c |}    " ///
		%8.6g CorrS_0 - CorrS_1 "   " %8.6g P_CORR "   "
	display "{text}{hline 12}{c BT}{hline 69}{c BT}{hline 23}"
	
	foreach result in X_0_mean X_0_var Y_0_mean Y_0_var XYcor_0 ///
		X_1_mean X_1_var Y_1_mean Y_1_var XYcor_1 ///
		CorrS_0 CorrS_1 ///
		N_0 N_1 S_0 W_0_mean S_1 W_1_mean VW_0 VW_1 ///
		VS_0 VS_1 P SP P_0 P_1 R_0_fisher R_1_fisher R_fisher ///
		P_CORR_0 P_CORR_1 P_CORR	 {
			return scalar `result' = `result'
		}
	
	//scalar list
	
	
end
